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Abstract 

Dynamical mechanisms that can stabiUze the coexistence or diversity in biology are generally 

of fundamental interest. In contrast to many two-strategy evolutionary games, games with three 

' strategies and cyclic dominance like the rock-paper-scissors game (RPS) stabilize coexistence 

Ph . and thus preserve biodiversity in this system. In the limit of infinite populations, resembling the 

Q I traditional picture of evolutionary game theory, replicator equations predict the existence of a 

■ fixed point in the interior of the phase space. But in finite populations, strategy frequencies will 

run out of the fixed point because of stochastic fluctuations, and strategies can even go extinct. 

For three different processes and for zero-sum and non-zero-sum RPS as well, we present results 

of extensive simulations for the mean extinction time (MET), depending on the number of agents 

^ ■ N, and we introduce two analytical approaches for the derivation of the MET. 

I> 
(N 

^ : 1 Introduction 

en 

^ ■ There are several examples of cyclic dominance in biology: Males of the califomian lizards Uta 

O ! stansburiana are known to inherit three different mating strategies which cyclically dominate 

each other Il22ll24ll . Another example is given by three strains of Escherichia coli |[T4l . in tropical 
marine ecosystems f?] or in high arctic vertebrates [10]. Cyclic dominance is also important 
^ \ in some theoretical models like the susceptible-infected-recovered-susceptible (SIRS) model in 

^ \ epidemiology 121, in cyclic extensions of the Lotka-Volterra model OTllSll or the Public Goods 

Game lfT3l[TTIl . Here, cyclic dominance is a way to preserve the coexistence of strategies (what is 
often called 'biodiversity' in a biological context) in the models. 

A straightforward model system for cyclic dominance is the rock-paper-scissors game (RPS) 
well known from schoolyards. It contains the three strategies 'rock', 'paper', and 'scissors' with 
a simple domination rule: Paper wraps rock, scissors cut paper, rock crushes scissors. The payoff 
a dominating strategy gets from the dominated one is set to 1 for normalization, the payoff for a 
tie is 0, and a dominated strategy gets a payoff —s, where we assume —s < 0. Hence the payoff 
matrix of this game reads 

/ -s 1 \ 
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For the standaid choice s = 1 we have a zero-sum RPS, for all other values the game is non 
zero-sum. One can simply intuit the impact of s Q : For huge s it is important for a player to 
avoid loosing, so it is successful playing the same strategy as the majority; hence an equilibrium 
that includes all three strategies is unstable. On the other hand, for small s it is more important to 
win occasionally so that the mixed equilibrium becomes stable. Experiments indicate s < 1 for 
the lizard example [36] and s > 1 for E. coli [14J- 

For a long time, the resulting dynamics have been described by a meanfield approximation in 
the limit of infinite populations. The traditional way to describe such evolutionary games was the 
standard replicator equation, an equation of motion for the density of an arbitrary strategy a [12], 

x„ = Va:a(vr"(x) - (7r(x))), (2) 

where is a constant prefactor which can be absorbed in the time scale by a constant rescaling of 
time. 7r°(x) = xiPq i + X2Pa,2 + X'^Pa^z and (vr(x)) are the payoff of strategy a and the average 
payoff of the whole population, respectively. The adjusted replicator equation |[23l 

K(x) - (7r(f))), (3) 



r + (vr(f)) 



has been used less frequently. The prefactor l/(r + (7r(x))) can be absorbed in the time scale by a 
dynamical rescaling of time for symmetric conflicts like the RPS, preserving stability properties. 
For asymmetric conflicts, like Dawkins Battle of the Sexes [7], this is not possible, as the average 
payoffs in each population do not coincide; hence the adjusted replicator equation can lead to 
qualitatively modified dynamics Il26l l6l. The derivation of the replicator equations and Fokker- 
Planck-equations (comprising the first-order corrections) for the intrinsic noise are commonly 
taken from the master equation of the underlying discrete stochastic process E7ll28ll26l[3Tl|29l . 

Both replicator equations predict the existence of a fixed point at (x/j, xp, xs) = (1) |) |)- 
The FP is neutrally stable for s = 1, asymptotically stable for s < 1, and unstable for s > 1. In 
the case of finite populations of size N, the population can drive out of the fixed point because 
of stochastic fluctuations, and after some time one of the strategies will go extinct. Once this 
has happened, a second strategy will die out soon because of the lack of cyclic dominance, and 
the third strategy will therefore survive forever. Several efforts have been made to overbear this 
shortcoming especially of the zero-sum RPS, for example spatial discretization of populations (for 
a review see |,25 |), mobility [21], the introduction of intelligent update rules (best response 131). 
the introduction of the parameter s as mentioned above 13, or the computation of a critical system 
size above which coexistence of strategies is likely f26 |, but it is still an open question how long it 
takes on average until the first strategy has gone extinct. In this paper we investigate the scaling of 
the mean time to the extinction of the first strategy (mean extinction time, MET), depending on the 
system size N and the parameter s, for well mixed populations and three evolutionary processes 
(with two of them having the standard and accordingly the adjusted replicator equation as limits 
N — )• oo), and present two analytical approaches which give theoretical insight in the scaling of 
the MET. 



2 Evolutionary processes 

Contrary to replicator equations describing the dynamics of relative abundance densities, real pop- 
ulations are finite (and discrete), appropriate modeling thus bases on discrete stochastic processes 
(of birth and death). In the classical Moran process [ 16] an individual a is chosen with probabil- 
ity proportional to its fitness. It reproduces, and the offspring replaces another randomly chosen 
individual b. In the frequency-dependent Moran process ifTSi which extends the classical Moran 
process |[T6l by considering coevolution, the fitness is not static but depends on the frequencies of 
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the strategies. For better comparison with the processes mentioned beneath, in each time step we 
choose an individual a at random which reproduces with probabihty 

6''^^\n n,n)- ^ ^ " ^ + ^M^) 

{na,np,n^) - 2 i - uj + u {tt) ' ^ ^ 

and the offspring replaces another randomly chosen individual b. Greek lettes a, /3, 7 can assume 
each of the strategies R, S, P, respectively. Here, 00 is the selection strength, the rij are the number 
of agents playing strategy i, -Ka = {'niPa,i+n2Pa,2+'nzPa,'i) / (-^—1) the payoff for an individual 
of strategy a, and (vr) the average payoff of the whole population. In the limit N — )• 00, the Moran 
process leads to the so-called adjusted replicator equation |26|. Note that a factor ^ is introduced 
here for a better comparability with the two processes mentioned beneath. 

The Moran process is a well-established stochastic process for evolutionary birth-death (for 
growing population sizes, see e.g. fSSl) dynamics with overlapping generations and therefore 
(in its frequency-dependent extension) serves as a standard model of evolutionary game theory. 
However, the Moran process requires perfect global information about the whole population, an 
assumption that can be unrealistic and undesirable. Because of that, two local processes have 
been mentioned. Again we choose an individual a at random for reproduction. Another randomly 
chosen individual h changes its strategy to the strategy of a with probability 
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■^lu {na,np,n^) = - H (5) 



for the local update |[26l and 

,l3{b)^a{a) , „ N \ (fi\ 

ct^F (n., np, n,) - ^ ^ ^_^(^^(^,_^^,,,) (6) 

for the Fermi process f34l, respectively, where Aurnax is the maximum of the possible payoff 
difference. In the limit N ^ 00, these processes lead to the standard replicator equation (local 
update) [26], and a nonlinear replicator equation (Fermi process) [30] with similar properties, 
respectively. The common approach for the derivation of the equations of motion and the first- 
order corrections for demographic stochasticity are to derive a Fokker-Planck equation from the 
master equation for the respective stochastic process I f27ll28ll26l l6ll29l. 

Here we focus on the mean extinction time for the Moran and Local Update processes. For 
each process, the probability of increasing the population strength of the strategy a by 1 in a single 
time step and decreasing the population strength of strategy /3 by 1, is given by 

Tprocess{na,ni3,n^ nQ, + 1, n/3 - 1, n^) = ^'^4>^^cessil^a,ni3,n^). (7) 

This quantity is known as hopping rate. Note that the sum over all hopping rates is 7^ 1 because 
reactions are possible that do not lead to changes in strategy frequencies. The time scale is chosen 
so that one reaction takes place every unit time step, including empty steps with no strategy 
change. 



3 Mean extinction times 

Compared to the time scale a mutant occurs - re-introducing a strategy - the time scale charac- 
terizing the survival properties of a genetic strain is defined by the mean extinction time. For two 
strategies, and fixed N, one has a onedimensional Markov process for which a closed expres- 
sion allows to proceed analytically as has been demonstrated by Antal and Scheuring UJ. For 
higher-dimensional cases unfortunately exact general solutions of the mean extinction time (or a 
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Figure 1 : Semi-logarithmic plots of the simulation data for the METs for the local update process 
(a), the Moran process (b) and the Fermi process (c). The METs were divided by N"^. Each point 
is an average over the extinction times of 1000 simulations. Red squares: s = 1.0, tu = 0.50, black 
rhombi: s = 1.0, to = 0.05, orange upwards triangles: s = 1.2, to = 0.50, yellow downwards 
triangles: s = 1.2, tu = 0.05, green open circles: s = 0.8, to = 0.05, blue filled circles: s = 0.8, 
CO = 0.50. For the Moran process cu = 0.45 has been used instead of u; = 0.50 to keep the transition 
probabilities in the interval [0, 1]. 

mean-first-passage time |[T9l ) problems are not known; in this case the situation is furthermore 
hampered by the location-dependent dynamics within a simplex boundary. So in most cases one 
has to rely on systematic numerical investigations. We have carried out extensive simulations to 
quantify the mean times until one of the three strategies has gone extinct. We have analyzed the 
mean extinction times for all three described processes depending on the number of agents N and 
the parameters s and lo. In general, we observe the following properties (see Fig.[Tl): For all three 
processes the MET becomes independent of the selection strength w for s = 1.0. Further on, the 
MET for s = 1.0 is identical for all three processes and proportional to N'^, 

(text)(s = l)~(0.54±0.02)iV2. (8) 

As expected, for s < 1 the MET is larger and for s > 1 smaller than for s = 1, but for weak 
selection the difference is small. In a logarithmic plot of the MET divided by N'^ one can easily 
see that {text{s = 1.0)) oc N"^, but also {text{s < 1)) oc Ar2e/(^,c^)JV ^jjj^ f{s,uj) > and 
independent of N. For the Moran process the stabilizing effect for s < 1 is larger than for the 
other processes. 

On the first view, for s > 1 the MET seems to be {text{s > 1)) oc N'^e^^^"''^^'^. But 
such a proportionality would imply a disappearance of the MET for large N. This is unrealistic 
because even the shortest way would require steps. For that a dependence {text{s > 1)) = 
ciN'^ + C2N'^e~^^^'^'>^ is more likely (with /(s, u) having the same properties as for the s < 1 
case), as it is predicted by one of our analytical approaches. 

In summary, it can be said that the MET switches from polynomial to exponential scaling at 
s = 1. This agrees with the drift reversal picture from Il26l l5l. So the (global) Moran process and 
the (local) local update and Fermi process behave similar, but the impact of the parameters s and 
uj is much stronger for the Moran process. 

4 Analytical approaches 

An exact analytical solution for two- or higher dimensional Markov chains is impossible, and a di- 
rect computation of the MET of such a three-strategy evolutionary game is not feasible. Although 
some efforts have been made, for example the work of Reichenbach et al. |f20l. wich have ana- 
lyzed mean extinction properties in an urn model of a three-strategy game, employing the usual 
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approach of applying van Kampen's linear noise approximation and deriving a Fokker-Planck 
equation [35, 27, 28, 26]. 

Adapting this approach to the zero-sum RPS for the local update and the Moran process, we 
find a predicted MET proportional to N"^ and independent of u, which is in accordance with the 
simulation data. Unf ornately this approach has some shortcomings: Although the overall scaling 
is correct, the numerical value (prefactor) of the predicted MET is too great by a factor of 5, and 
it is not possible to use this approach for the non zero-sum RPS or in the Fermi process. For that 
we have developed two analytical approaches which do better. We will present both approaches 
for the local update. Following the same schemes for the other two processes is, although possible, 
a bit more sophisticated because not all integrals can be done analytically here in general. As we 
have seen from the numerical investigations, the s / 1 payoffs override differences between 
the underlying processes, so that it is warranted to concentrate on one analytically more feasible 
process in the remainder. 



4.1 First approach - Expected changes in the distance to the fixed 
point 

To compute - with help of approximations - the MET for the general case, we need an appropriate 
coordinate system. For the standard replicator equation and s = 1, 

H = -XRXpxs (9) 

is a constant of motion which assumes the value H = —1/27 in the fixed point and = on the 
edge of the phase space, the simplex 5*3. Here, xr, xp and xs are the frequencies of the strategies 
R, P and S, respectively. For s < I, H is a Lyapunov function of the replicator equations with 
H < 0, and so the inner fixed point ist asymptotically stable [5|. For s > 1 the fixed point is 
unstable, and the attractor of the system moves to the edge of the simplex. Via a transformation 
of the fixed point to the origin of the phase space, and by inserting xr + xp + xs = 1, which 
guarantees the conservation of the total density, we find for H: 

H = -{x + l/3)iy + l/3){l/3-x-y) (10) 

For large absolute values of H, the curves with H = const, resemble slightly deformed circles, 
but for — they approach the triangle form of the simplex. In the following we will use H as 
an observable to measure the effective distance to the fixed point. But obviously we need a second 
coordinate to describe a point in the phase space, so we will use a system with the variables x and 
H by eliminating the y-coordinate. This gives us two solutions for y{x,H), an upper branch 
above the fixed point, and a corresponding lower branch below. We need both because it is not 
possible to describe all points of the phase space (x, y) by only one solution (see figure [2(b)] ). So 
we have 

_ -3x - _ ^4 + 108F + 122; + 2,2AHx - 27x2 - 54x3 + Slx^ 

6(1 + 3x) ^^^^ 

and 

-3x - 9x2 _^ 74 + i08i? + 12x + 324i?x - 27x2 - 54x3 + §1^4 

■ ^ ^ 

Here the indices 'u' and 'd' stand for 'up' and 'down', respectively, because by re-inserting of x 
and H in the first (second) equation one will get only values of y that lie above (under) a separator. 
In two points, Xmin and Xmax, both solutions are identical so that the curve H = const, is closed. 
We can compute these from the condition yu = yd and solving for x. 
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(a) (b) 



Figure 2: (a) A plot of curves for which H = const, holds, for different values of H. Different 
colours on each curve mark the different solution branches ?/„ und yd. The black line connects the 
values Xrnin ^ud Xrnax, respectively, (b) Plots of the simulation data for the MET from the local 
update (symbols) for different pairs of s and cu depending on N and the corresponding results from 
the analytical approach presented in section l4~n (lines . with a not filled version of the corresponding 
symbol at the end). Red circles: s = 1.0, cu = 0.5, orange squares: s = 1.2, to = 0.5, yellow triangles: 
s = 1.2, cu = 0.05, green rhombi: s = 0.8, to = 0.05. 



i cos Q arg (^-54H + QV?>yjH{27H + 1) - l)^ 
-\/3sm Qarg {-UH + Q\f?,yjH{27H + 1) - l)^ + 1^ 



(13) 



and 



cos Q arg (^-54i/ + 6V^VH{27H + 1) - l)^ 
+\/3 sin Q arg (^-MH + QV?,yjH{27H + 1) - ^ + 1 



(14) 



and a third (formal) solution which here is of no use because for typical values of H it produces 
values outside the simplex. arg(2;) is the argument of the complex number z. 

Now we compute the expectation value of the change in H for a single time step for the local 
update process. For time t, the system may be in an arbitrary allowed point (x, H). For that we 
will need the hopping rates in the (x, H) coordinates. For example, for the change 5\ 



AT' ' 



' N 



the hopping rate reads after inserting y = Uu, 



T 



with 



1 



1 



N 



-9x2 _^ 3^ _^ J _^ 2) 
216(3x + 1) 



(27^x2 + 9(^ - 2)x + (2s + 1)/V' - 6) (15) 



/ = \/(3x + 1) ((3x + 1)(2 - 3x)2 + 108i7). 



(16) 



Note that we will use all variables as if they were continous, although in fact they are not, and 
we use T(5i,52) as an abbreviation for T(xi,X2,X3 — )• xi + 6i,X2 + 82, x-^ + ^3) due to the 
conservation of total density (where 5i und 82 can take all values of {— 0, as long as 
they are not identical). 
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Figure 3: (a) A plot of the expectation values of the changes in H, 5H{H), multiplied with A^^, for 
the local update process. The dashed lines show this quantity for = 798, the continous lines for 
N = 396. red : s = 1.0, uo = 0.5, orange: s = 1.2, uj = 0.5, yellow: s = 1.2, uj = 0.05, green: 
s = 0.8, CO = 0.05, blue: s = 0.8, co = 0.5. (b) Additional illustration of the meaning of 5H{H) 
(s = 0.8, CO = 0.50). If the observable h has a certain value Ht at a time t, H is expected to make a 
step in the direction of the arrow belonging to this specific value of H. The lengths of the arrows are 
over-emphasized by a factor of 40000. 



Further on we will need the change in H that emerges from a change of the (x, y)-values in 
the (x, i/) -coordinates as well. To derive this, we write H{x + Si,y + 62) of the changed values 
of X and y, subduct the initial value H{x, y), and transform the result into the (x, i/)-coordinates. 
For y = Vu and the same change as in the example above we get 

/I 1 \ i-21Nx'^ -%N + 2)x + Ng-Q) . ^ ^ , 

with 



g = yj{2,x + 1) (27x3 - 27x2 + 108i7 + 4) (18) 

By multiplying all hopping rates T {51,82) with the associated changes in H, /S. {51,82), and 
summing over all terms, we get a relatively simple equation for the expectation value of the 
change in i7 in a single time step, 

5H{x,H) = Y,T{5i,52)^{5i,52) (19) 
H{9 + 27x + (1 + 27H)N{-1 + s)^^ + 27iV(-l + s)x^^^) 

with 5ilu (x, K) = 5Hd (x, H) =: 5H (x, H). Averaging over x by integrating over x from Xmin 
to Xmax and normalizing by the length of whole interval we derive 

(6a2 + (66 - 3)a + 66^ -36 + 2) N{s - 1)^ + 18 



5H{H) 



6A^2 

(20) 



(a - 6)iV 

with a = Xmin, b = Xmax- The simplest way to approximate the mean extinction time for a 
system starting in the fixed point of the replicator equations is then to define a map 

H,+i = Hi + 5H{Hi) (21) 
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and to iterate it until H(j;^ is reached. The MET is then simply the number of necessary steps. 



As one can see in the right hand side of Figure 2(a) this result is in good agreement with the 
simulation data, for the dependence of s and w, as well as for the numerical values. For s = 1 the 
MET is predicted to be proportional to N'^ and independent of uj because 6H{H) is independent 
of a; for s = 1. It is in good agreement with the simulations that the MET for s 7^ 1 and weak 
selection differs only slightly from the MET for s = 1, while for strong selection and s > 1 it is 
clearly smaller than for s = 1. 

Unfornately, for s < 1 this approach works only for small selection and not to great A'^, 
because for s < 1 non-positive expectation values of the changes in H are possible, so if we start 
the mentioned iterated map in the fixed point and repeat the iteration again and again, we would 
arrive at a value of H where the expected change is at some time, and the predicted MET would 
be infinite. The part of the interval where the expectation value is negative is bigger for larger uj 
than for smaller uj, and it grows with increasing N, so in this cases this approach is not sufficient. 

4.2 Second approach - Fokker-Planck equation 

For the local update process, a Fokker-Planck equation (EPF) for the time evolution of the proba- 
bility density is given by (summation over double indices) 



dtPix, t) = -d^[aiix)P{x, t)] + -di,[Bijix)Pix, t)], (22) 



where the coefficients are 

ai 

and 



6xiT{x ^x + 6x) (23) 



Sx 



Bij{x) = 5xi5xjT{x — )■ X + 5x), (24) 
which gives the following results 



Sx 



{2xp + XR){l + 3xR)i; 

OtR = — (25) 

2,N ^ ' 

(1 + 3xp)(xp + 2xr)V> 
ap = — (26) 

2 + 'ixR- 

BrR = ^^^2 (27) 

2 + 3xp -9x?, 



8 



Transforming this into polar coordinates xr = r cos 0, xp = r sin ^, we find for the coefficients 
of theFPE 

r^(3;- cos(o) + 1)(2 cos(a) + 2(.s + 1) siii(o) + r{s - l)(Hiii(2o) + 2)) 
= 6iV 

(30) 

rV'lSr sin(^) + l)(-2(s + 1) cos(0) - 2ssin(^) + r(s - l)(sin(2^) + 2)) 

6iV 

(31) 

-9r^cos^((/.) + 3rcos((/)) + 2 
Brr = ^^^2 (32) 

(3rcos((?!)) + l)(3rsin((^) + 1) 
Srp = BpR = ^ (33) 

-9r^sin^((/>) + 3rsin((/.) + 2 
= ' (34) 

where = ui/AiTmax, and hence for the FPE itself 

p(OAi) = 6iV(. - 1)^ sm{24>y + N (l2r^ - l) (. - 1)^; - 9 (o,o,o) ,,,, 

37V2 ^ 
_^ /rV'(-3(s + 2) cos ((?!)) + (7s + 2) cos (3^) - 2(4s + (2s + 7) cos (2^) + 5) sin(^)) 



12A^ 

(s + l)V'(sin(20) + 2) cos(0) + 3 cos(3(^) - sin((/>) + 3 sin(30) cos(2(/)) \ (o,i,o) 

m 12iVV + Qjsf2j.2 J 

rV((2s + 1) cos(^) + (2s + 7) cos(30) - (s + 2) sin(^) + (7s + 2) sin(3^)) 



i2Ar 

r^{s - l)V'(cos((/)) sin(0) + 1) r{-N{s - + N{s + 1) cos(20)^ - 18) 



N 

cos{(f)) - cos(3^) + sin(^) + sin(3^) cos(0) sm{(p) + 1 \ (i,o,o) 
8iV2 ^ 9N^r J 

_ (cos(0) - sin(0))(9cos(0) sm{(j))r + 3r + cos(0) + sin(0)) ^(^^ ^^^g) 

9r cos(^) - 9r cos(3(?!)) + 9r sin(0) + 4 sin(2^) + 9r sin(3(?!)) + 8 p(o,2,o) 

— 36r^ + 3 cos(^)r + 9 cos(30)r + 3 sin(0)r — 9 sin(3^)r — 4 sin(20) + 8 ^2 0) 
^ 72iV2 " ' 

withP(^''^-) :=|^|,|.P(r,0,t). 

If we now assume the probability density to be independent of the angle (which is a 

good approximation at least for small r), all terms which contain derivatives with respect 
to (f) will drop out. By averaging over cf) afterwards by integrating from to 2% over and 
dividing the result by 2% we get a probability density which is only independent of r and 
t. Hence the FPE reads 

^ r(N(12r'^-l)(s-lU-9)^, , 
dtP(r, t) = — — ^— ^-P(r, t) (36) 



^3(iV(6r--l)(.-l)^-18)r- + 2 



2 

2 ) 



=: a{r)P{r, t) - v{r)drP{r, t) + D{r)d;P{r, t). (37) 
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-2. X 10"^ : 
-4. X 10"^ : 
-6. X 10"^ : 



Figure 4: Left: A plot of the convection velocity t>(r) for the local update with = 3960. As one 
can see, v{r) is ~ for s = 1. This corresponds to the neutral stability of the fixed point in the 
limit — 7- oo for s = 1. For s < 1 f (r) is negative, corresponding to the asymptotical stable FP, 
and for s > 1 v{r) is positive corresponding to the unstable FP. For r — t- v(r) diverges for all 
parameter values to — oo, which results from the simplifications in combination with the existence of 
a FP at r = 0. Right: Semilogarithmic plot of simulation data of the MET in the local update process, 
divided by A^^ (symbols) for several pairs of s and tu, and the corresponding values predicted by the 
approach in section l42l (lines, with a not filled version of the corresponding symbol at the end). Red 
circles: s = 1.0, tu = 0.5, orange squares: s = 1.2, tu = 0.5, yellow rhombi: s = 1.2, to = 0.05, 
green upward triangles: s = 0.8, oj = 0.05, blue downward triangles: s = 0.8, tu = 0.5 



Let us now approximate the diffusion by its value at r = 0, this gives us 



D 



1 



(38) 



r can theoretically vary between and 2/3 for some values of 0, but most extinction 
events will happen around those point of the edge of the simplex that are closest to the 
inner fixed point. These points have the distance r > ^ from the inner FP. Because of 
that we search for the solution of the one-dimensional random walk descripted by eq. [37] 
in the interval [0, L] = [0, |]. Next, let us approximate the convection velocity t>(r) by its 
value for mediate r, e.g. i" = ^ = ^, and we get 



-0.00462963AfV)-0.166667 
Af2 ' 



]_ 

0.00462963jV^-0. 166667 

JV2 ■ 



0.8 
1.0 
1.2 



(39) 



(We cannot approximate it by its value at r = as we have done for the diffusion because 
f (0) is — oo corresponding to the FP.) For a last simplification, we neglect the term that 
includes the probability density itself. Hence, the FPE reads 



dtP{r, t) = -vdrP{r, t) + DdlP{r, t). 



(40) 



The computation of the MET for such an equation with one reflecting border at r = 
and one absorbing border at r = L is a standard problem, for example, see lfT9l . One can 
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find the MET by introducing a time integrated probability density, 

Po{r) = / P{r,t)dt. (41) 
Jo 

With this, the time-integrated FPE reads: 

= (42) 

or ox^ 

with the boundary conditions Po(^) = Oatr = L and J(r) = j{r,t)dt = vP{r,t)- 
DdrP{r,t)dt = 1 at r = 0. There, J(r) is the time-integrated probability current. The 
constraint J(0) = 1 corresponds to the injection of an unit probability current at t = 0, 
this means a start of the system in the origin. The solution of the differential equation is 
given by 

pJr) = - (1 - e"(-^)/^) (43) 

V 

Then the MET is simply 

(t) = r P,{T)dT =^-^[l- e-'^^/^), (44) 
Jo V 



what leads to 



or in general 



Af2('~0.00154321Ari/;+0.183191eO-0i38889^'/' -0.16666?) 



S = 0.8 



(0.00462963Afi/'+0. 166667)2 ' 

0.594885Ar2^ s = (45) 

-0.oi38889iV^^2('gO.0i38889iVV<(0.00154321Ar,/.-0.166667)+0.18319l) 



(0. 166667-O.OO462963Af-0)2 



1.2 



72e-4^(--i)^-5Ar2 (^e*^(«-i)'^+l(5Ar(s - 1)^ - 108) + 72e) 
^ (36 - 5N{s - • ^"^^^ 

The proportional dependence of the MET on N and the other parameters is very good, 
as one can see in Fig. |4l This approach predicts the MET to be proportional to A^^ 
and independent of the selection strength a; for s = 1, while for s < 1 it forecasts the 
MET to be substantially proportional to A^2gKii/>Af^ s > 1 in good approximation 

proportional to ciN'^ + 026"^'^^. 



5 Summary 

Evolutionary games with three cyclically dominating strategies can stabilize the coexis- 
tence of strategies superior compared to games with only two strategies which include 
a dominating strategy. In finite populations, it is still likely that one of these strategies 
goes extinct, though on average it takes a time which is considerably longer than in two 
strategy games. We have analyzed the mean times to the extinction of one of the three 
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strategies numerically in the RPS-game for three different evolutionary processes, the lo- 
cal update, the frequency-dependent Moran process and the Fermi process. These mean 
extinction times (MET) have the same fundamental dependencies of the number of agents 
A^, the parameter s and the selection strength cu. 

For the zero-sum RPS game (s = 1), the MET is proportional to N"^ and independent 
of u! for all three processes; even the proportionality constants are the same. For s < 1, we 
find a stabilization of biodiversity, as the MET grows exponentially with A^. In contrast to 
that, for s > 1, there is a destabilization of biodiversity, and the MET now grows smaller 
than in the zero-sum RPS (where it grows oc N'^), namely, a part of the proportionality 
factor of the zero-sum RPS now decays exponentially with A^. In both non zero-sum 
cases, the coefficient that is multiplied with the N in the exponent, grows with both u 
and (s — 1). 

Solving such two-dimensional Markov chains is not possible, but we have developed 
two analytical approaches for the approximation of the dependencies of the MET of N 
and the parameters s and to. The first approach can be used for s > 1 only, but serves 
for all three processes. In contrast, the second approach gives a good approximation for 
all parameter values, but it can be used completely analytical only for the local update 
process, as for the other two processes not all integrals can be done analytical. 

In summary, the fixation time - as a long-time property of the process - shows a 
crossover from exponential to polynomial scaling with the population size, being fully 
consistent with the critical population size delved from the "drift reversal" picture which 
is based on the short-time dependence of the average drift in phase space - repelling 
versus attracting towards the fixed point. Thus both approaches to describe stochastic 
stability in finite populations here lead to the same conclusions: coexistence is preserved 
(stabilized) for a positive-sum cyclic game, whereas negative-sum games as well as small 
populations destabilize coexistence. 
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